Maximum size of drops levitated by an air cushion 
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Liquid drops can be kept from touching a plane solid surface by a gas stream entering from 
underneath, as it is observed for water drops on a heated plate, kept aloft by a stream of water 
vapor. We investigate the limit of small flow rates, for which the size of the gap between the drop 
and the substrate becomes very small. Above a critical drop radius no stationary drops can exist, 
below the critical radius two solutions coexist. However, only the solution with the smaller gap 
width is stable, the other is unstable. We compare to experimental data and use boundary integral 
simulations to show that unstable drops develop a gas "chimney" which breaks the drop in its 
middle. 



I. INTRODUCTION 

Drops levitated on an air cushion have numerous ap- 
plications, and have generated interest for a long time. 
For example, in lens manufacture drops of molten glass 
can be prevented from contact with a solid substrate [l| . 
This is achieved by levitating the glass above a porous 
mold, through which an air stream is forced. A second 
example is the so-called 'Lcidcnfrost' drop 0, a drop of 
liquid on a plate hot enough to create a film of vapor be- 
tween the drop and the plate d 0, H Q . Since the drop 
is thermally isolated insulated by the vapor film, it can 
persist for minutes || ■ Finally, a thin air film is believed 
to play a crucial role for the "non-coalescence" of a liquid 
drop bouncing off another liquid surface @, H, ■ 

The question we will address in this paper is whether 
for a given set of parameters, in particular the radius of 
the drop as it "rests" on the substrate, a stationary so- 
lution exists and whether it is stable. Apart from lense 
manufacture [l|, this question is important for the ma- 
nipulation of corrosive substances [Kj or the frictionless 
displacement of drops [||. Of particular interest is the 
maximum drop size that can be sustained, and the limit 
of very small flow rates. The drop continues to levitate in 
this limit since the gap between the liquid and the sub- 
strate becomes very small, so the lubrication pressure 
produced by the viscosity of the gas becomes significant. 
This enables us to employ asymptotic methods, making 
use of the disparity of scales between the gap size and 
that of the drop. 

Experimentally, it is observed that the stability limit 
is reached when the radius equals at least a few capillary 
lengths i c = y/j] (pg). This natural length scale for our 
system is determined by the surface tension 7, density p 
of the liquid, and acceleration of gravity g. At a few cap- 
illary lengths, the drop is flattened to a pancake shape. 
Biance et al. d EJ observed a critical radius 



4.0 ±0.2, 



(1) 



where r max is defined in Fig. [TJ Beyond this radius, 
"chimneys" appeared, i.e. bubbles of air trapped below 
the curved and concave surface of the drop that rise ow- 
ing to buoyancy, and eventually burst through the center 
of the drop. This suggests that the critical radius is re- 
lated to the Rayleigh- Taylor instability of a heavy fluid 
the drop) layered above a light fluid (the gas layer). In 
l~fl ] this idea is used to estimate r max /£ c w 3.83. 
While this is close to the experimental value, the argu- 
ment ignores the gas flow responsible for the levitation 
force. This flow was taken into account by Duchemin et 
al. who calculated the static shape of a drop lev- 
itated above a curved porous mould, using a combina- 
tion of numerics and asymptotic arguments. For large 
drop volume, they found no physical solutions, while for 
smaller drops multiple solutions were calculated numeri- 
cally. Duchemin et al. Ill] did not present a formal sta- 
bility analysis of these solutions, but suggested that the 
limit of stability is related to the appearance of large os- 
cillations on the underside of the drop. 

A large number of studies of Leidenfrost drops have fo- 
cused on the appearance of self-sustained oscillations of 
the drop @, H O, El EI El EE E3- These oscillations 
can sometimes lead to a morphological bifurcation of the 
drop which takes the shape of a star [1 El El EE E3 • 
Similar star-shaped drops have been reported in drops 
vertically vibrated on non-sticky surfaces, and the shape 
is generally attributed to a parametric instability [18l ] . 
Our original question was whether oscillations could per- 
haps be explained even in the limit of viscous drops, 
which we focus on in this paper. This is not the case, 
since both our asymptotic results and simulations of the 
complete dynamics (i.e. beyond linear stability analysis) 
show that once unstable, a drop breaks up owing to the 
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formation of a chimney. 

We treat both the liquid drop and the surrounding 
gas in the incrtialcss (Stokes) limit. For the asymptotic 
analysis, we also require the drop to be much more vis- 
cous than the gas. The main effect of this assumption is 
that there is hardly any flow inside the drop, so it can 
be treated as being in hydrostatic equilibrium at any in- 
stant in time. We also prescribe the rate at which gas 
is injected into the underside of the drop, thus ignoring 
the possible interplay between drop dynamics and vapor 
production in the Leidenfrost problem. 

Our analysis is similar in spirit to the earlier paper of 
Duchemin et al. [l[ , but we only address the simpler case 
of a flat substrate. As a result, we are able to perform 
all the calculations analytically (up to a few universal 
constants, which have to be computed numerically). Our 
solution curves are in qualitative agreement with those 
for a curved substrate [l| , but now imply a full analytical 
description. In addition, we are the first to perform a 
systematic stability analysis of the stationary states. We 
find the maximum stable radius 




inner neck region 



FIG. 1: Definitions and sketch of the matching regions. 



flux Q as the relevant quantity, which can be calculated 
by integrating the gas flux entering from underneath up 
to the neck position r„ . Let us also introduce a slightly 
different dimcnsionalization of the flow rate 



4.35 



(2) 



where f goes to zero in the limit of vanishing gas flow. 

For typical experimental flow rates we find that f sa 
0.4, consistent with the experimental result (p}. At the 
end of the paper we discuss how our analysis relates to 
the stability argument of [H, [ll[ , based on the Rayleigh- 
Taylor instability. 



II. PROBLEM FORMULATION 



A. Geometry and dimensionless parameters 
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which will appear naturally in the analysis. 

Finally, another parameter is the viscosity ratio be- 
tween liquid and gas 



A 



_ ??drop 
— : 
Vgas 



(6) 



but which will be considered asymptotically large for 
most of this paper. Throughout, lengths will be ex- 
pressed in £ Cl velocities in 7 /rj gaai and stresses in "f/i c - 



We consider axisymmetric drops of liquid, levitated 
above a flat surface by gas flowing into the underside 
of the drop, cf. Fig. [1] We set out to find the shape of 
stationary drops and their stability, as a function of the 
gas flow rate and the drop volume. The size of the drop 
is expressed by the Bond number 



(3) 



where V is the volume of the liquid drop, and R = 
(47rl73) ^ the unperturbed radius. The dimensionless 
gas flow rate supporting the drop reads 



QVg 



Pel 



(4) 



where Q is the volume of gas that escapes through the 
narrow neck region (see Fig.[T]) per unit of time, and 7y gas 
is the viscosity of the gas. Our analysis will identify the 



B. Structure of the problem 

The problem we set ourselves is to solve the inertialess, 
axisymmetric fluid flow equations, with a prescribed in- 
flux of gas into the underside of the drop. Most of our 
analytical work assumes in addition that the drop is much 
more viscous than the gas. The structure of the expected 
solution is shown in Fig. [TJ The gas pressure below the 
drop has to be sufficiently large in order to support the 
weight of the drop. In the limit of small dimensionless gas 
flux r, the gap between the drop and the substrate must 
therefore be small in order to generate enough pressure. 
The underside of the drop inflates to a gas pocket, whose 
width is of similar size as the drop itself. The narrow 
gap is formed in a small neck region only, where a large 
curvature assures that the gas pressure can be sustained 
by corresponding surface tension forces. Apart from this 
viscous neck region the gas pressure is constant, both in 
the gas pocket as well as to the exterior of the neck. 

This leads to the following asymptotic structure of the 
problem, characterized by the matching between three 
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different regions. In the limit of small flux all viscous ef- 
fects become localized in a small neck region, situated at 
a radius r — r n from the center. In this region, there ex- 
ists a balance between viscous and surface tension forces. 
In addition, the slope of the gap profile h(r) turns out to 
be small in this region, so lubrication theory permits 
to reduce the flow equations to an ordinary differential 
equation for h(r). We will call this the inner solution or 
neck region. 

To close the problem, boundary conditions are needed. 
These are provided by two outer regions on either side 
of the neck, denoted by '— ' (the gas pocket toward the 
center of the drop), and '+' (the outside of the drop). 
Both regions are controlled by a balance of gravity and 
surface tension alone. First, we solve the equations in 
each of the regions individually. Second, we require that 
both the slope and the curvature of the profile match 
smoothly at the boundaries between two regions. This 
leads to a set of equations which determines stationary 
drop solutions uniquely. Solutions exist only below a 
certain critical neck position r c , in which case we find two 
branches, one with a small gap width (the lower branch), 
and an upper branch with a larger gap width. 

Our stability analysis of the two branches is based on 
the observation that the relevant dynamic variable is the 
position r„ of the neck, which can shift easily. The max- 
imum stable neck radius does not coincide with r c , but 
is significantly smaller, located on the lower branch. 

III. INNER SOLUTION: NECK REGION 

A. Lubrication approximation 

We consider incompressible, axisymmetric flow in the 
gas layer, so that mass conservation gives 

rh + {rhu)' = rv(r). (7) 

Here u is the depth averaged horizontal velocity of the 
air layer, while v{r) is the rate at which air volume is 
injected per unit area below the drop. The main focus of 
the paper will be on stationary states and their stability. 
Stationary drop profiles are found by taking h = 0, and 
integrating ([7|) to 

rhu=^, (8) 

where T(r) = 2ir dr'r'v(r') is the flux in the lubri- 
cation layer. In the case where the injection source is 
localized at r = 0, the flux T is simply constant. 

To get a closed equation for h(r) in the neck region, 
we solve for u. As our results will confirm, the neck 
region is shallow, ti <C 1, meaning that we can use the 
lubrication approximation [l9| to analyze the flow, see 
Fig. [1] Owing to the large viscosity ratio between the 



drop and the surrounding gas, the liquid drop acts as a 
no-slip boundary, and the flow in the gas layer is well 
approximated by 




Since the Reynolds number is very small in typical ex- 
periments, we use the Stokes approximation [19j to relate 
this velocity to the pressure. As there is almost no flow 
inside the drop, the liquid will be at hydrostatic equilib- 
rium, pnquid = Po — z - Furthermore, the pressure jump 
across the interface equals the curvature times the surface 
tension, so one obtains the lubrication pressure inside the 
gas layer as 



p = p -h-h". (10) 

In what follows we will show that the width of the neck 
region scales as T 1 / 5 and thus is asymptotically small in 
the limit of vanishing flux. We are therefore permitted to 
neglect the axisymmetric contribution to the curvature in 
the neck region. Using the horizontal component of the 
Stokes equation, p' = d 2 u/dz 2 , we find 

u = 12h 2 (ti + ti") . (11) 

Now © , (|11|) provide a closed equation for the stationary 
interface profile h(r): 



h 3 (ti + ti") = ^^. (12) 
7rr 

The right hand side of (fT2|) represents the viscous stress 
in the flow, and will only become important when h is 
small, i.e. in a small neck region around r n , where we 
may set r — r n . This gives 



ti (ti + ti") = x , (13) 

with 



x , (») 

7TTY, 

A crucial observation is that there is no need to know 
the precise form of how the gas is injected, but one only 
requires the flux across the neck. This of course provides 
a great simplification for the Leidenfrost problem, where 
evaporation rates are related in a complicated way to the 
temperature profile inside the drop. 
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B. Similarity solution for neck region 

As gravity is unimportant in the thin neck region, (|13p 
can be further simplified to 

h 3 h"' = x- (15) 

Since we are interested in the limit of small flux, we look 
for similarity solutions 

h(r)= X a H(0, where £ = (16) 

giving 

H 3 H"' = 1, with 4a - 3/3 = 1. (17) 

In the limit £ — > oo, the solutions have to match onto a 
sessile drop of constant curvature. Since 

h" = X a ~ 2P H", (18) 

one needs that a — 2/3 = for the curvature to remain 
finite as x ~~ * 0- Together with (fl7|) this fixes a = 2/5 
and /3 = 1/5, and hence we have 

Hr) = x^B {^f) • (19) 

The form of the similarity function will be determined 
from the matching below. The fact that a > (3 justifies 
the assumptions made so far. First, we find that h' <C 
1 in the limit x ~ > 0, justifying the use of lubrication 
theory. Similarly, h! <C h'" , so that both gravity and the 
axisymmctric curvature can indeed be neglected in the 
neck region. 

The asymptotic behavior of (jTTJ) is quadratic for both 
£, -> ±oo, 

H+ = ±K + ? + S+t for e^cx) (20) 
#_ = ]^K-i 2 + S-i for i -> -oo. (21) 

Physically, the values of the asymptotic curvatures K± 
set the pressure in the corresponding outer regions. 

Since (JTTJ) is of third order, solutions can be specified 
by three independent parameters, one of which can be 
absorbed into a shift of £. Therefore, the two asymptotic 
curvatures K± uniquely determine the solution. As a 
consequence, the slopes S± are dependent variables. To 
perform the matching, we require the function 




FIG. 2: Solid line: the function / relates the slope S- to cur- 
vature K'_ of the inner solution (Eq. (|23p with K+ = 2.17, 
corresponding to r n — ro). Dotted lines: the function g pro- 
vides the matching condition between inner region and gas 
pocket region (Eq. (|51|l with x = 10~ 7 ). The three curves 
correspond to values r n = 3.55 (below critical), r n = 3.62 
(critical), and r n = 3.65 (above critical). 



whose existence is assured by the above argument. Since 
(|17p is invariant under the transformation h — > h/a and 
x — > x/a 4 / 3 one must have 

f(K.,K + )=K 1 + /5 f(^-,l\. (23) 

This function is computed numerically and is plotted in 
Fig. [21 We show below that stationary solutions cor- 
respond to the intersection of / with another function 
<7, shown in the same figure. It can be seen that the 
matching breaks down at a critical neck radius r n , be- 
yond which stationary solutions cease to exist. 



IV. OUTER SOLUTIONS 

Having seen that viscous effects are localized in the 
neck region, the rest of the drop is at static equilibrium. 
Hence, the pressure is constant both in the gas pocket 
between the drop and the substrate, as well as to the 
exterior of the neck. These pressures are not equal, how- 
ever, since one requires a pressure difference to drive the 
flow across the neck. In Fig. [T] we therefore distinguish 
two outer regions, denoted by + and - respectively. Since 
P± = Pa^h — K, the outer solutions can be obtained from 

K + h = c±, (24) 



S-=-f(K_,K+), 



(22) 



where k is the curvature of the interface. The constants 
c± determine the pressure difference across the neck 



FIG. 3: The outer solution of the "drop" region corresponds 
to a perfectly non-wetting sessile drop. The size of the drop, 
characterized by r n , sets the curvature K+ for the inner so- 
lution. 



FIG. 4: The outer solution fixes the value of K+ as a function 
of the neck radius (which controls the drop volume). The 
maximal radius ro = 3.38317 ■ ■ ■ gives K + = 2.17 •■ ■. 



Ap = p_ -p + = c + -c_, (25) 
and will follow from the matching. 

A. Outer "drop" region: + 

Below we will find that the profile of the "drop" re- 
gion requires dh/dr — > as h — > 0, in order to match 
to the neck smoothly. This corresponds to a perfectly 
non- wetting sessile drop (Fig. [3j. When matching the 
curvature we also require d 2 h/dr 2 = K + as h — ► 0. Ow- 
ing to the vanishing slope near h = we are allowed to 
write k = d 2 h/dr 2 in (|24[) . Hence, one finds c+ = K + . 

To deal with the overhang of the sessile drop, it is 
convenient to solve the profile in terms of the arclcngth 
s along the interface. We define 9 as the angle with the 
horizontal and rewrite (1241) as 



~dl 

dh 

ds 
dr 

ds 



sm I 



- h + K, 



sm ( 



— = cos#, 



with boundary conditions 

6(0) 
h(0) 
r(0) 
0(s t ) 
r(s t ) 






r„ 

7T 

0, 



(26) 
(27) 
(28) 



(29) 
(30) 
(31) 
(32) 
(33) 



where St is the value of the arclcngth at the top. Two of 
these five boundary conditions serve as the definitions of 
r n and st, so that the remaining three boundary condi- 
tions fix the solution uniquely. The equations have been 
solved numerically. 

Each value of K + thus gives a solution with a differ- 
ent r n , some of which are shown in Fig. [3] The nu- 
merically obtained relation between K + and r n is de- 
picted in Fig. |4j For the maximal neck radius r n = ?*o = 
3.38317 ■ • •, introduced below, one finds K + = 2.17 ■ • •. 

The value of r„ effectively sets the volume of the drop. 
Namely, the weight of the sessile drop is carried by the 
pressure exerted by the substrate on the contact area 7rr„ . 
The difference between the liquid and the gas pressures 
at h = is simply K +1 so we find 



K+irr 2 , = 2nV+ 



V+ 



1 



(34) 



where V+ is defined as the real volume divided by 27r, i.e. 



dh ixr . 



(35) 



1 

Note that to obtain the real liquid volume, one has to 
subtract the volume V- of the gas pocket. However, V- 
goes to zero in the limit of vanishing gas flow, as shown 
below. 



B. Outer "gas pocket" region: - 

In the "gas pocket" region, the profile h(r) is no longer 
multivalued and we can express the curvature as 



h" 



(1 + ^2)3/2 r (l + /j'2)l/2 



(36) 
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The solution is then specified by (|24|) with boundary con- 
ditions 



h'(0) = 
h(r n ) = 0. 



(37) 
(38) 



In Appendix [Al we show that the solution can be written 
as an expansion 



h(r) 



Jp{r) - Jo(r n ) 
Jo(r n ) 



+ 0(c 3 _) 



(39) 



where Jo(r) is a Bessel function of the first kind. Using 
furthermore that the curvature has to match the curva- 
ture of the inner solution K-, and thus K- = h"(r n ), we 
can further simplify to 



h(r) = K_ 



M?) - Jo{r n ) 



0{Ki 



(40) 



We see that the thickness scale of the gas pocket is set 
by the value of K_. In the limit of vanishing flux we 
expect this thickness to tend to zero, making K_ a small 
parameter. To find solution branches, it is crucial to go 
beyond linear order and to find the term of order K^_ in 
(|4U|) . The only quantity that is needed to perform the 
matching to (|22[) . coming from the inner solution, is the 
slope h'_(r n ). This calculation is done in Appendix |A"1 

At this point we can already infer an upper bound on 
the possible values of r n . Figure [5] shows the outer gas 
pocket solution (with normalized amplitude) for various 
values of r n . The outer solutions are defined on the do- 
main where h(r) > 0, hence the maximum possible neck 
radius is achieved when Jq(t) has its first minimum at 
the maximal radius ro = 3.8317 - • • (vertical line). The 
corresponding solution is drawn with a heavy line. 

A second remark is that Jq vanishes at r w 1.852 ■ ■ -, 
so K- must become zero at this radius. At even smaller 
radii Jq turns negative, which yields negative values of 
K-. However, the inner solution cannot reach large h for 
negative K_ , which means that the matching procedure 
described here does not work. Dealing with this problem 
requires an additional matching region between the inner 
and outer (gas pocket) solution, the introduction of which 
lies beyond the scope of this paper. We will simply stay 
away from r n « 1.852 and instead focus on radii close to 
the maximal value ro ~ 3.8317, as detailed below. 



MATCHING THE ASYMPTOTIC REGIONS 



FIG. 5: Outer solutions for the gas pocket region (amplitudes 
normalized to unity). Thin solid lines correspond to r n = 1, 2; 
the dashed line, corresponding to r„ = 3, illustrates that h 
would have to become negative to realize a neck radius larger 
than ro . The heavy line shows the maximum possible r n = ro , 
corresponding to the first minimum of Jo(r). 



hout± = ^±k(r-r„) 2 +/4k(r-r„) (41) 
h in ± = \K ± {r-r n ) 2 + X 1/b S±{r-r n ). (42) 
Therefore the matching conditions become 



X 



K± 



h'l\r„ 
h'± \r„ 



(43) 
(44) 



The conditions on the curvature were already taken 
into account when computing the outer profiles from 
Typical values for K + are of order unity, while the 
slope requires h' + \ rn = as \ ~~ * 0. This is why for the 
first outer solution we considered a perfectly non- wetting 
drop. 

The — conditions are more subtle. The thickness of the 
gas pocket goes to zero asymptotically so that both h!'_\r n 
and h'_\ Tn will be small. In this case the selection of the 
solution explicitly requires the slope condition, which we 
express as 



S- = 



K_ h'_ 



-g(K-,r n ;x)- 



X 1/5 h'L\r n 

Together with (|22|) this closes the matching problem 



(45) 



A. Matching conditions 

We can now match the asymptotic regions by express- 
ing (|20[) , (f2"Tj) in their original variables and expanding the 
outer solutions around r = r.„ , 



f{K-,K + )=g{K-,r n ; X )- 



(46) 



This equation indeed contains the three matching re- 
gions: K + implicitly depends on r„ through the + outer 
solution, / is determined by the inner solution, while g 
follows from the — outer solution. 
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B. Bifurcation: critical radius r c 



For a given value of the flux Xi we have reduced the 
problem to finding the intersections of the functions / 
and g. This is sketched in Fig. [21 showing / and g for 
X = 10 -7 and several values of r n . Depending on the 
value of r n , there can be two intersections, one intersec- 
tion (when the curves are tangent), or no intersection. 
Each intersection corresponds to a stationary drop solu- 
tion. This can be translated into a bifurcation diagram 
showing K_ vs r n (Fig. [7]). For small radii there are two 
branches of solutions, corresponding to the two intersec- 
tions, which merge at r c . No stationary drop solutions 
exist for r n > r c . 

We analyze the bifurcation in the limit of vanishing 
flux, x "~ * 0. We will show that 



r c = r + O(x 



,2/15 



(47) 



so that the critical neck radius r c approaches the maximal 
radius r$ in the limit of vanishing flux. To analyze the 
vicinity of the critical point we introduce 



r = r - r n . 



(48) 



At the same time we will find that K_ oc x - This 
means that as the limit of x going to zero is reached, 
K- = and r n = rg, which implies K + = 2.17 according 
to Fig.|H These two data fix the solution of (|T7"|) uniquely, 
and lead to the asymptotic profile shown in Fig. [5] From 
its minimum, one finds that 



FIG. 6: The inner solution H(£) obtained from numerical 
integration of (T7}, with K+ = 2.17 and K- = 0. It will 
follow that these values correspond to the critical solution. 
The minimum value H n « 0.931 determines the thickness of 
the neck (HTfl). 



For details we refer to Appendix \Kl where we show that 
g 2 = 1.486 

The matching condition f = g (cf. (I46]) ) is now re- 
duced to a horizontal line intersecting a cubic function: 



Y 



Solving for f , one finds 



■h = —( f -92K 2 _). 



(53) 



0.931 x 



2/5 



(49) 



in agreement with the scaling found by [l[. 

We now analyze the first correction to the solution as 
X increases, but in the limit where x,-ff_,f -C 1. This 
can be done by considering the corresponding limit of 
the functions / and g, cf. Fig. [51 Namely, the function / 
approaches a constant, which is found numerically to be 

/~/o = 1.12.-.. (50) 

on the other hand, the asymptotic form of g becomes 



9- X 



- 1/5 K_ (f - g 2 K 2 



(51) 



The first term of (|5ip is found by expanding (|40|) for r n 
close to r : 



Jo(ro)(r-r n ) 



h'L\ 



-f + 0(K 



0{K 2 _) 



(52) 



where we used the property Jq(^o) = 0. We need to keep 
the K 2 term as it can become of the same order as f . 



r(K-, X ) = 



A/5 f 
C 1 Jo 



(54) 



which has been plotted for different values of x in Fig. [7] 
Thus for a given value of r n one finds two solution 
branches, which end at the critical value 



1/3 



V 



2/15 + q 



(x 4/15 ) 



(55) 



as claimed before. Note that the smallness of the power 
2/15 makes r c non-negligible. For typical experimental 
values of the flux the critical point is thus substantially 
shifted with respect to the asymptotic value tq. 

Plugging this back into (f54|) . one finds the value of K- 
at the critical point: = 0.72x 1/15 . But it follows 

from (H that h = h(Q) m K_(l - J (r ))/Jd(r ), and 
thus the maximum gap width is to leading order 



2.52 X 



1/15 



(56) 



This concludes the analysis of the stationary solutions, 
which are described by (|B"4")) . At a given flow rate the 
critical neck radius is given by (|55|) . which approaches 
the maximal value tq in the limit x ~~ * 0- 
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FIG. 7: The bifurcation diagram (K-,f), derived from 
(|54[) . Curves correspond to different values of the flux, \ = 
10 -5 , 10 -7 , 10 -10 , revealing the weak dependence on \- The 
dashed lines represent perturbations Sr,8K- = —Sf/c which 
are tangent to the solution curve. They represent marginal 
perturbations, separating stable from unstable solutions. 



VI. STABILITY 

We now turn to the important question of which part 
of the solution branches shown in Fig. [7] are stable. Es- 
sentially, we find that the lower branch is linearly stable, 
while the upper is linearly unstable. Surprisingly, how- 
ever, the marginal point is not exactly at the maximum 
radius, but slightly before. 



of Ap. Had the pressure stayed constant, F would be 
larger than the weight of the drop leading to the for- 
mation of a chimney and thus to instability. Similarly, 
pressures smaller than the marginal condition leads to a 
stable situation, giving the stability criterion 



Ap' < 0. 



(59) 



In the limit of small x, the pressure difference (|25l) is 
simply the difference of the curvatures: 



Ap = K + - K- , 
so that stability requires 



K'-K 1 



2Ap 



< 0. 



(60) 



(61) 



The derivative K' + can be read off from Fig. 2J and is 
negative. Clearly, this has a stabilizing effect. The sign 
of K'_ can be inferred from the bifurcation diagram. The 
lower branch has a stabilizing contribution, while the 
upper branch is be destabilizing. The location of the 
marginal point, however, depends on the numerical val- 
ues of the three terms. 

Taking the derivative of we find 



K' 



2V+ 



2K, 



(62) 



A. Stability limit 



Moreover, for vanishing flux A'_ <C K + , hence we may 
replace Ap ~ K+ , giving the stability criterion 



Once more we make use of the fact that there is no flow 
inside the drop, so that the drop shape adjusts quasi- 
statically to variations in the neck region. We therefore 
consider infinitesimal variations in the neck position, 6r n 
and assess the corresponding change in levitation force 
SF. Since the pressure difference Ap = p_ — p + across 
the neck acts for r < r n , this force reads 

F = Apnrl. (57) 

A marginal perturbation Sr n occurs whenever the re- 
sulting levitation force is unchanged, 5F = F'5r n = 0, so 
that it still equilibrates the weight of the drop. Hence, 
we find the marginal condition 

(58) 

where the prime denotes the derivative with respect to 
r n . In order to produce the same levitation force, an in- 
crease in r n thus has to be compensated by a decrease 



K'_ > (63) 

' n 

Near the maximal radius r n w r the criterion for stabil- 
ity becomes 

2V 

K'_ > c- 1 = —± = 0.92 ■ • • . (64) 
r a 

Indeed, the upper branch with K'_ < is unstable, but 
the marginal point is not at the maximum radius K'_ = 0, 
but slightly before. This is indicated in Fig. [7] by the 
dashed lines, which each have a slope 0.92. 

B. Linear stability analysis 

We now include dynamics into the stability analysis. 
We note first that an infinitesimal variation of the neck 
position, 5r n = —5r, also induces a variation of the cur- 
vature, SK_, and of the flux, S\. These three parame- 
ters are related through mass conservation of the liquid 
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and the gas. The analysis is closed by a third equation 
coming from matching the dynamic inner region to the 
hydrostatic outer regions. 

The volume of the liquid can be calculated by the vol- 
ume enclosed by the sessile drop solution, V+, minus the 
volume of the gas pocket V-, i.e. 



liquid = V+(r n )-V-(r n ,K-). 



(65) 



This is exact up to asymptotically small corrections due 
to the inner region. The volume V+ is (numerically) de- 
termined by the value of r n , while V_ can be computed 
analytically using (f40|) . 



drrh-{r) 



-K r 2 

2 n 



J (r„) - 2J 1 {r n )/r r 
Jo(r n ) - J\{r n )/r n 



(66) 



The expression simplifies at ro, because of the property 
Ji(fo) = 0. Since the liquid volume is strictly conserved, 
liquid = 0, one finds near r 



5K_ = -- 



Sr 



(67) 



where the constant c has been defined by ([64)) . Relation 
([6"7j) expresses the fact that when r n increases, increasing 
V + , the volume of the gas pocket has to increase by a 
similar amount to keep the liquid volume constant. This 
is achieved by an increase of K_. 

Mass conservation of the gas is described by continuity 
(jTj) , which can be integrated to 



addition the slope of such a tangent curve must be — c _1 
according to (|67p . this uniquely fixes a point on any of the 
lines at constant \. The critical tangent curve is drawn 
dashed in Fig. [71 Below this point, on the lower branch, 
solutions are stable, above they are unstable. 

Formally, the growth rate of perturbations is computed 
by writing 



Sf = a8r. 



(71) 



Now using (fBTj) . ([70)1 . and the first variation of ([S"4"[) one 
finds 



(72) 



24 


fdf\ 


-l 


df 


cro 


[dx) 









The partial derivatives are to be evaluated from ([54)1 . 

This indeed gives the same stability boundary as (1631) . 
which was based on a global force balance (note that 
dr/dK- = — (A'i) -1 ). The maximum stable radius r s is 
found by the condition a = 0, yielding 



X 1/5 /o 
K 2 



2g 2 K_ = c 



(73) 



as an equation for K_ at the stability boundary. This 
value of K— is inconsistent with the asymptotic estimate 
K_ sa x 1 ^ 15 considered so far, indicating that the point 
where the solution exchanges stability is at a distance 
slightly larger from the critical point r c . This means 
that A"_ is smaller than expected (further down the lower 
branch, cf. Fig. [7]). Thus the second term on the right of 
(|73p is small compared to the other two, and we obtain 



rhu = r(r) 



d_ 

ot 



dr r h. 



(68) 



The second term on the right hand side can be identified 
as the rate of change of gas pocket volume V- , which we 
will write as —Vl_6f. This change absorbs part of the 
injected air, decreasing the flux passing across the neck. 
Considering the radius somewhere inside the neck region, 
r ?s r n , the equation can be simplified to (using (jlip and 
h' < 1): 



h 3 ti" = X + 5X, 
where the variation of the flux reads 



(69) 



to\ 1/10 



K- = [~J X 1/w - (74) 
If we evaluate (|54[) in the same limit, we finally obtain 



(foc) 1/2 X 1/W . 



(75) 



Thus for vanishing flux the maximum radius of stable 
solutions approaches r , but with an even smaller power 
than r c . This scaling implies that r s < r c < ro, as seen 
in Fig. [3 



VII. NUMERICAL TESTS 



(70) 



The matching condition (|54j) closes the dynamical sys- 
tem, taking into account the dependencies (f6"T|) and ([70)) . 
The marginal case 8\ — corresponds to a curve tan- 
gent to any of the lines r(K-) shown in Fig. [7J Since in 



A. Non-linear dynamical behavior 

We begin with a simulation of the full axisymmetric 
Stokes problem, using a boundary integral method [2fj| . 
which has the advantage that it tracks the interface with 
high precision. The idea is to regard the interface as a 



continuous distribution of point forces, which point in 
the direction of the normal and whose strength is pro- 
portional to the mean curvature. Since for Stokes flow 
one knows the Green function giving the velocity field 
resulting from a point force, one can write the velocity 
anywhere in space as an integral over the free surface. In 
an axisymmctric situation, the angle integral can be per- 
formed, so the remaining integration is one-dimensional. 

External flow sources can simply be added; in the 
present case we take the gas flow as a point source of 
strength Q situated at the origin on the solid plate which 
bounds the flow. For this a simple exact solution is avail- 
able [H]]. Likewise for the Green function one must take 
into account the presence of a no-slip wall. This is pos- 
sible using the method of images [22||, and the resulting 
boundary integral formulation has been applied success- 
fully to the motion of drops relative to a wall [23| . If, as 
in our case, the viscosity of the drop is different from that 
of the surrounding, one must account for the stress mis- 
match across the interface. This can be done at the cost 
of introducing another integral over the velocity on the 
interface into the equation, which turns the equation for 
the velocity field into an integral equation. After solving 
this equation for a given interface shape, the thus com- 
puted velocity field can be used to advance the interface. 

We follow closely an earlier implementation of the 
boundary integral method, used for example the coales- 
cence of two drops inside another fluid [24J. The only 
significant difference is that the free space Green func- 
tion has been replaced by those for half space, bounded 
by a wall. We tested the code by comparing to an ex- 
act solution of a sphere moving perpendicular to a wall 
[2l| . This is realized in the limit of a very small drop, 
or of very large drop viscosity, so that there is hardly 
any deformation. The agreement was good, but signifi- 
cant deviations occurred when the gap between the wall 
and the drop was smaller than 5 % of the drop radius. At 
present, we do not know the origin of this numerical prob- 
lem, which prohibits us from investigating the asymptotic 
limit of very small gap spacings. Instead, we report on 
simulations at moderate gap spacings, which show the 
nonlinear stages of chimney formation, not captured by 
our linear stability analysis. 

Figure [5] shows a viscous drop which is slightly smaller 
then the stability boundary. Starting from a configura- 
tion shown as the light curve, it relaxes toward a station- 
ary stable state (heavy line). For a Bond number which 
is just slightly larger, the same initial condition leads to 
a rising gas bubble in center of the drop, see Fig. A 
thin film forms between the rising gas bubble and the top 
of the drop, which drains slowly. As seen in Fig. the 
neck radius is r« ~ 2.5, giving \ = 0.015. On the basis 
of our asymptotic theory ([7S]) . a rough estimate of the 
stability boundary gives r s w 3.2. 
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FIG. 8: Boundary integral simulation of a drop with parame- 
ters T = 0.02, Bo=4.2, and A = 100. The drop relaxes toward 
a stable state, which is drawn as the heavy line. 
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FIG. 9: The same as Fig. [8] but with a slightly larger Bo = 
4.4. The air bubble under the center of the drop lifts up to 
form a chimney. The time interval between the profiles is At 
= 3000, in units of £ c %as/7- 



B. Lubrication approximation 

To test the bifurcation scenario in more detail we re- 
sort to direct numerical simulation of the lubrication 
equation. Due to the overhang of the drop, we sepa- 
rate the upper part of the drop and the lower part of 
the drop at the maximum radius, r max , defined by the 
point \h'\ = oo. The upper part is solved as described in 
Sec. IIVI A. and for the lower part of the drop we use 



K (1 + h«)3/2 + r (l + h*) 1 / 2 1 ' 

X = h 3 ( K ' + ti). (77) 
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This describes both the inner and outer regions in the 
lower part of the interface, while we have conveniently 
taken the rate of injection T(r)/r to be constant for all 
r. Boundary conditions for this 3rd order equation are 



h'(0) 







^(^"max ) ^patch? 



(78) 
(79) 
(80) 



where K pa t c h is the curvature at the point where the upper 
and lower solutions are patched. A one-parameter family 
of solutions is obtained through variation of the upper 
part of the drop. It was shown in [l[ that this procedure 
provides drop solutions that arc quantitatively accurate. 

The numerically obtained drop profiles are conve- 
niently characterized by the position of the neck, r n , and 
the gap below the center of the drop, h . Numerical re- 
sults for the solution branches are shown as solid lines in 
Fig. [TO] for two values of the flux. \ = 10 -4 is a typi- 
cal experimental value encountered for Leidenfrost drops, 
while \ = 10~ 7 illustrates the convergence toward the 
asymptotic limit. As predicted, there is a critical radius 
beyond which no stationary solutions exist. The asymp- 
totic predictions shown in Fig. [7] have been translated to 
the dashed lines of Fig. [TO] These are obtained from ([54"]) . 
where K_ was computed from ho using (|40[) . Good quan- 
titative agreement is achieved for small enough values of 
the flux. 

Finally, we determined the critical radius r c for a range 
of values of the flux \. Figure [TT] shows how the numer- 
ical values (dots) indeed approach the asymptotic pre- 
diction (solid line) in the limit of vanishing flux. Due 
to the very small powers x 2 ^ 15 : the convergence towards 
ro = 3.8317 • • • is extremely slow (horizontal line). As a 
consequence, the correction with respect to this asymp- 
totic value will be significant for typical experimental val- 
ues of the flux. 




FIG. 10: Bifurcation diagram ho vs r„ for \ — 10 -4 and 10 -7 . 
Smaller \ yield larger radii. Solid lines were obtained from 
numerical solution of the lubrication equation (|77|l . Dashed 
lines correspond to asymptotic theory (|54ll . 
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VIII. DISCUSSION 

Owing to the smallncss of the neck region (|4"9"]l . we can 
make the simplification that the pressure inside the gas 
pocket below the drop is constant (Fig. Q~| . This pressure 
is larger than the atmospheric pressure and provides the 
force required to levitate the drop. Matching the pressure 
difference across the neck with the viscous flow then pro- 
vides the bifurcation diagram of Fig. [7] yielding a critical 
neck radius r c . In the limit of vanishing flux, the criti- 
cal radius approaches ro = 3.8317- ■ •. This value arises 
because it is the first minimum of the function charac- 
terizing the shape of the gas pocket, which is the Besscl 
function Jo(r). For larger r n , the gas pocket shape would 
need to become negative, which is of course not allowed. 

Experimentally, the size of the drop is measured by 
looking at the drop from above. This measurement pro- 
vides the maximum radius r max rather than the neck ra- 
dius, cf. Fig.[TJ For large puddles the difference between 



FIG. 11: Critical radius r c as a function of the flux \- 
The numerical values obtained from the lubrication equa- 
tion (dots) indeed approach the theoretical prediction (|55[) 
(solid line). The dashed line indicates the asymptotic value 
r = 3.8317 



r max and r„ approaches V2 — arcccosh-s/2 « 0.53. For 
drop sizes relevant here we confirmed numerically that 
''max — T n ~ 0.52. Combined with (|55p . we thus find 

W, st « 4.35 - 1.02 X 1/10 (81) 

for the boundary of stability, expressed in terms of the 
capillary length. Typical experimental values of x can 
be extracted using (|49p . and typical experiments yield 
h n fa 100/i to, obtained from diffraction data This 
gives x ~ 10 -4 , and thus r maXjSt ~ 3.95, to be compared 
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to reported experimental values of 4.0 ± 0.2 0, [TT|. A 
similar estimate of \ is obtained from considering the la- 
tent heat of evaporation [f|. Furthermore, our boundary 
integral simulations show that the nonlinear dynamics 
for larger drops lead to the formation of chimneys, as ob- 
served experimentally. We are therefore confident that 
the analysis in terms of Stokes flow provides an accurate 
description of this instability. 

Let us now return to the argument put forward in 
[S Ej|, relating chimney formation to the Raylcigh- 
Taylor instability. The later occurs when a layer of fluid 
is suspended above another fluid of lower density, so that 
the system tends to destabilize due to buoyancy forces. 
Surface tension opposes this effect, so that the instabil- 
ity occurs at long wavelengths only. Biance et al. [ll| 
propose that levitated drops remain stable as long as ax- 
isymmetric perturbations that fit inside the drop are sta- 
ble with respect to this buoyancy driven instability. 

For an infinitely extended liquid film, one finds that 
Jo(fcr) are axisymmetric eigenmodes, with the stability 
criterion k > 1. While the Bessel function does not have 
a well-defined period, the maximum drop size was esti- 
mated in [ll1 | by the first minimum of the mode with 
k = 1, occurring at tq. In hindsight, our results justify 
this choice of taking the minimum of Jo(^) as the sta- 
bility boundary, provided that it is identified with the 
neck radius, rather than with r max . With this connec- 
tion, our results reduce to the Rayleigh- Taylor argument 
in the limit of vanishing gas flow, showing that the bal- 
ance between buoyancy and surface tension provides the 
right mechanism. The effect of the gas flow is to slightly 
reduce the range of stable solutions (|8Tj) . 



APPENDIX A: GAS POCKET SOLUTION 

In this appendix we expand the gas pocket solution 
for small amplitudes and compute the constant <?2- We 
consider the equation 



h" 



(1 + ^2)3/2 r (i +h ny/2 
with boundary conditions 



h = 



(Al) 



y 



y 



+ y = o, 



( 1+y /2)3/2 r (l+y'2)l/2 

with boundary conditions 

2/(0) = 

y(r n ) = -A, 
where A = c_ We expand in A, 

v( r ) = Ayi(r) + A 3 y 3 (r) + O (A 5 ) 
This yields a hierarchy of equations 



(A4) 



(A5) 
(A6) 



(A7) 



y't + ^ + yi 

r 

ii , 2/3 , 

% + — + yi 

r 



0. 



(A8) 
(A9) 



with boundary conditions 



yi{r n ) 

2/3(0) 

2/3 {r n ) 



(A10) 
(All) 
(A12) 
(A13) 



The first equation gives yi(r) = —Jo(r)/Jo(r n ), which 
can be inserted into the right hand side of the equation 
for y 3 (r). 

To compute the constant 32, we require the ratio 
y' ( r o) / y" ( r o) ■ I n terms of the expansion 



2/M 
y"(ro) 



A 2 y' 3 (r ) + O (A 4 



(A14) 



where we used the properties yi(ro)' = and ^"(ro) = 
— 2/1 ( r o) = 1- Comparing to (|51[) . we simply find 



h'(0) 
h{r n ) 

This is equivalent to solving 



(A2) 
(A3) 



92 = 2/3 M- 



(A15) 



We obtained this value numerically by solving the ODE 
for 2/3, for which we numerically obtained 92 = 1-486 ■ • •. 
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